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INTRODUCTION 


This  is  the  second  yearly  technical  report  on  a  three  year 
effort  to  study  physical  processes  of  relevance  to  the  mass 
spectrometric  measurement  of  stratospheric  ions.  The  effort 
involves  the  development  of  a  Monte  Carlo  model  of  the  free  ;:et 
expansion  occurring  within  the  mass  spectrometer  including  the 
effects  of  agglomeration  onto,  and  fragmentation  of,  ionic 
clusters . 

The  attempt  to  carry  out  in  situ  mass  spectrometry  in  the 
stratosphere  is  complicated  by  changes  that  may  occur  in  the  gas 
stream  as  it  expands  after  passage  through  the  orifice.  Both 
positive  and  negative  ions  exist  in  the  stratosphere  with  clus¬ 
tered  polar  molecules  surrounding  the  ion  core.  As  these  ion 
clusters  are  carried  along  in  the  expanding  gas  stream,  the 
falling  temperature  will  tend  to  favor  the  formation  of  larger 
clusters.  The  charge-dipole  interaction  is  characterized  by 
large  cross  sections,  so  agglomeration  of  polar  molecules  may 
change  the  cluster  size  distribution  that  the  quadrupole  sees 
from  the  distribution  that  exists  in  the  undisturbed  stratos¬ 
phere.  Conversely,  the  measured  cluster  size  distribution  may 
be  driven  towards  smaller  clusters  via  fragmentation.  As  the 
ionic  clusters  are  selectively  accelerated  by  the  electric  field 
within  the  mass  spectrometer,  high  energy  collisions  with  neutrals 
may  break  apart  the  clusters. 

The  present  effort  involves  a  Monte  Carlo  simulation  of 
these  processes,  so  that  a  model  can  be  used  to  relate  the 
measured  properties  to  those  existing  in  the  undisturbed  atmos¬ 
phere.  The  basic  elements  of  the  direct  simulation  Monte  Carlo 
method  were  described  in  the  previous  yearly  report,^-  and  they 

^Elgin,  J.  B.,  "Monte  Carlo  Calculations  of  Mass  Spectrometer 
Flow,"  Report  AFGL-TR-83-0057 ,  Air  Force  Geophysics  I.abontory , 
February  1983.  ADA  128069. 


will  not  be  repeated  here.  There  were  substantial  advancements 
made  in  the  past  year  which  are  described  in  detail  in  this 
report.  Section  2  describes  the  ability  to  treat  separa  .e  spa¬ 
tial  segments  of  the  flowfield  in  separate  computer  runs,  pro¬ 
viding  a  substantial  increase  in  computational  efficiency.  The 
formalisms  for  handling  the  simulation  of  internal  molecular 
energy  and  its  interchange  with  the  translational  mode  have  been 
generalized  in  terms  of  the  classic  Chi-Square  distribution, 
and  this  generalization  is  discussed  in  Section  3.  Means  of 
efficiently  sampling  from  a  Chi-Square  distribution  were  devel¬ 
oped,  and  are  discussed  in  Section  4. 

The  final  new  development  for  the  past  year  is  the  inclusion 
of  the  accelerating  effect  of  electric  fields  on  charged  species, 
and  this  item  is  discussed  in  Section  5.  Section  6  presents  a 
summary  of  the  code  status  and  Section  7  discusses  a  test  case 
calculating  the  diffusive  separation  of  a  binary  mixture  of  C02 
and  H2. 

2.  SPATIAL  SEGMENTATION  OF  SOLUTION  REGION 

The  code  was  generalized  in  the  past  year  so  that  it  has 
the  ability  to  compute  sequential  spatial  segments  of  the  solu¬ 
tion  starting  from  the  orifice.  (This  new  ability  is  merely  an 
option  and  in  no  way  affects  the  capability  of  handling  the 
entire  flow  field  at  once.)  The  use  of  this  option  offers  the 
potential  for  a  substantial  decrease  in  the  memory  and  computing 
time  required  to  solve  a  problem.  Details  of  the  spatial  seg¬ 
mentation  scheme  are  discussed  in  the  subsections  below. 

2.1  Justification  for  Segmentation 

The  first  question  that  must  be  addressed  is  whether  the 
segmentation  of  the  solution  is  physically  and  mathematically 


justified,  and  this  question  is  critically  related  to  the  ques¬ 
tion  of  boundary  conditions.  The  physical  laws  that  are  embodied 
in  the  Monte  Carlo  solution  procedure,  involving  molecular  trans¬ 
lations  and  collisions,  are  as  valid  for  a  portion  of  the  solution 
region  as  they  are  for  the  whole  region.  In  order  to  carry  our 
the  solution  in  just  a  subregion,  however,  it  is  necessary  that 
the  boundary  conditions  can  be  specified  a  priori  alone  >e 
boundaries  of  the  subregion  in  question. 

For  a  Monte  Carlo  flow  field  calculation,  the  boi  ary 
conditions  are  imposed  by  specifying  the  velocity  dist.  i  .tion 
function  for  incoming  molecules  along  all  boundaries,  and  then 
selecting  molecules  from  this  distribution  with  the  proper  frequecy 
and  introducing  them  into  the  simulation.  Usually  this  involves 
extending  the  boundaries  to  a  region  of  undisturbed  (known)  flow, 
or  to  a  region  where  molecular  backflow  into  the  solution  region 
is  insignificant. 

In  the  simulation  of  mass  spectrometer  flow,  the  upstream 
boundary  is  taken  to  correspond  to  one  dimensional  sonic  condi¬ 
tions  at  the  orifice,  except  that  the  mass  flow  is  reduced  by  an 
empirically  determined  discharge  coefficient.  The  solution 
region  is  extended  far  enough  downstream  and  to  the  side  so  that 
the  flow  of  molecules  into  the  solution  region  from  these  other 
boundaries  can  be  neglected.  (An  exception  to  this  is  the 
incursion  of  background  gas  into  the  jet  which  will  be  treated 
as  an  equilibrium  gas  at  the  side  boundaries.  This  feature  has 
not  yet  been  added  to  the  code,  but  it  does  not  affect  the  pre¬ 
sent  discussion  since  the  distribution  function  for  this  gas  will 
be  assumed  known  at  the  side  boundaries.)  Hence,  a  well  defined 
solution  can  be  carried  out  downstream  of  the  orifice  as  long  as 
the  downstream  boundary  is  far  enough  from  the  orifice  to  assure 
that  molecular  backflow  is  negligible.  Since  molecules  have  a 


8 


thermal  velocity  which  is  on  the  order  of  the  speed  of  sound, 
backflow  becomes  negligible  when  the  local  Mach  number  is  large 
compared  to  unity.  (It  is  interesting  to  contrast  this  to  a 
continuum  calculation,  which  requires  only  that  the  local  Mach 
number  be  greater  than  one,  perhaps  by  a  very  small  amount,  for 
there  to  be  no  upstream  influence.)  In  a  Maxwellian  gas,  the 
probability  of  an  individual  molecule  having  an  upstream  velocity 
direction  is  given  by  P,  where 

P  =  jerfc  (M  V772)  /  (1) 

and  M  and  y  denote  Mach  number  and  ratio  of  specific  heats, 
respectively.  For  sonic  flow  this  probability  is  on  the  order 

of  5%,  but  by  the  time  the  Mach  number  has  become  2.0  the  proba- 

-4 

bility  is  on  the  order  of  4x10  and,  for  all  practical  purposes, 
backflow  can  be  ignored. 

Hence,  the  imposition  of  the  downstream  boundary  condition 
neglecting  backflow  into  the  solution  region  is  justified  very 
shortly  after  the  orifice,  since  the  flow  rapidly  becomes  sub¬ 
stantially  supersonic.  As  long  as  this  is  the  case,  it  is 
perfectly  proper  to  solve  for  a  small  segment  of  the  solution 
region.  Once  that  solution  is  completed,  then  the  derived 
velocity  distribution  on  the  downstream  boundary  defines  the 
required  upstream  boundary  condition  for  the  next  solution  seg¬ 
ment.  This  information  is  automatically  written  to  a  file  which 
is,  in  turn,  automatically  read  as  input  for  the  next  segment. 
Segmentation  can  be  used  in  conjunction  with  the  separation  of 
major  and  minor  species,  if  desired. 

2.2  Advantages  of  Segmentation 

Although  the  physical  justification  for  segmenting  the 
solution  region  has  been  demonstrated,  there  remains  the  question 
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as  to  why  one  would  want  to  "turn  one  problem  into  two  or  more 
problems."  Frequently  the  subdivision  of  a  physical  problem 
into  multiple  subproblems  reduces  the  total  effort  involved,  and 
this  is  no  exception.  Specific  advantages  are  enumerated  below. 

2.2.1  Reduced  Storage 

When  the  problem  is  broken  up  into  segments,  the  total 
number  of  cells  and  molecules  required  in  any  one  segment  is 
less  than  required  for  the  larger  single  solution.  This  trans¬ 
lates  directly  into  a  decreased  requirement  for  computer  core. 

2.2.2  Decreased  Time  to  Achieve  Steady  Flow 

The  Monte  Carlo  technique  is  inherently  a  procedure  which 
solves  an  unsteady  flow  problem.  For  cases  such  as  the  present 
one  where  the  problem  of  interest  is  really  a  steady  state  flow, 
this  is  solved  by  letting  the  unsteady  solution  relax  to  a  steady 
state.  The  computation  time  required  to  achieve  a  steady  state 
can  be  regarded  as  "computational  overhead"  for  the  present 
problem,  since  useful  steady  state  sampling  of  the  solution 
cannot  begin  until  steady  state  has  been  achieved. 

The  time  required  to  establish  a  steady  flow  can  be  esti¬ 
mated  a  priori  as  the  length  of  the  solution  region  divided 
by  the  orifice  flow  velocity  (times  a  safety  factor) .  Hence, 
the  longer  the  initial  solution  segment,  the  greater  is  the 
period  of  unsteady  flow.  If  the  entire  solution  region  is  cal¬ 
culated  at  once,  then  the  program  may  well  spend  most  of  its 
effort  in  simply  achieving  steady  state.  If  the  solution  region 
is  segmented,  however,  then  the  first  segment  can  reach  steady 
state  substantially  faster  than  the  whole  solution  region  does. 
This  is  of  particular  significance  since  the  solution  near  the 


orifice  is  the  most  collision  dominated,  requiring  a  large  por¬ 
tion  of  the  computational  effort  needed  in  each  time  step.  When 
subsequent  segments  are  solved,  they  still  require  a  relatively 
long  time  to  achieve  steady  state  (though  the  time  is  somewhat 
diminished  since  the  solution  region  is  shorter) ,  but  the  compu¬ 
tational  effort  per  time  step  is  substantially  less. 


2.2.3  Increased  Time  Step  for  Latter  Segments 

The  time  step  in  a  Monte  Carlo  simulation  should  generally 
be  small  compared  to  the  mean  time  between  collisions  for  a  mole¬ 
cule,  since  the  processes  of  translation  and  collisions  are 
treated  separately  in  each  time  step.  If  the  entire  solution  is 
solved  for  at  once,  this  implies  that  the  entire  solution  is 
constrained  to  the  relatively  small  time  step  required  by  the 
collisional  region  near  the  orifice.  If  that  region  is  solved 
for  separately,  then  subsequent  segments  can  have  a  substantially 
larger  time  step  since  the  mean  time  between  collisions  is  larger 
in  subsequent  segments.  Hence,  even  though  the  latter  regions 
still  require  a  relatively  large  amount  of  simulation  time  to 
achieve  steady  state,  this  time  is  more  easily  accomplished 
since  the  allowable  time  step  is  much  larger. 


2.3  Segmentation  Summary 

These  considerations  strongly  suggest  that  the  first 
spatial  segment  should  be  made  within  a  few  diameters  of  the 
orifice,  since:  1)  The  flow  is  by  then  already  sufficiently 
supersonic  to  justify  the  neglect  of  backflowing  molecules,  and 
2)  The  number  of  time  steps  required  to  achieve  steady  state  is 
small,  which  is  particularly  important  in  this  collision  domin¬ 
ated  region  of  the  flow.  Note  that  the  question  is  completely 
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unrelated  to  whether  or  not  the  flow  is  in  equilibrium  (as  would 
be  the  case  if  an  initial  region  were  to  be  calculated  by  the 
method  of  characteristics) .  It  is  simply  a  matter  of  making 
the  calculation  of  the  flow  field  more  efficient  by  separating 
the  relatively  long  relaxation  time  which  is  required  by  the 
latter  portions  of  the  flow  field  from  the  collision  dominance 
which  is  characteristic  of  the  initial  portion.  Once  the  first 
segment  is  calculated,  then  subsequent  segments  benefit  computa¬ 
tionally  from  the  ability  to  take  substantially  larger  time  steps 
in  those  regions. 

3.  IMPROVEMENTS  IN  INTERNAL  ENERGY  MODELING 

3.1  Discussion 

In  the  previous  yearly  report,1  the  modeling  of  internal 
energy  of  molecules  and  the  interchange  between  the  internal  and 
translational  modes  was  discussed  at  length.  In  the  past  year, 
significant  progress  was  made  in  unifying  and  streamlining  the 
procedures  for  describing  molecular  energies ,  and  a  module  was 
developed  which  is  substantially  more  efficient  in  describing 
equilibrium  flow.  These  practical  improvements  are  fundamentally 
related  to  the  theoretical  result  that  all  of  these  cases  can  be 
related  to  sampling  from  a  Chi-Square  distribution.  The  theo¬ 
retical  developments  are  discussed  in  this  section,  and  means  of 
sampling  from  a  Chi-Square  distribution  are  discussed  in  the 
following  section. 

3.2  Initial  Internal  Energy  Levels 

The  classical  distribution  of  internal  energy  among  mole¬ 
cules  with  v  degrees  of  freedom  (Eq.  (86)  of  Ref.  1)  can  be 


properly  simulated  if  each  individual  molecule  is  assigned  an 
internal  energy,  Ej ,  given  by 

Ej  =  jRqTX  ,  (2) 

where  Rq  is  the  universal  gas  constant,  T  is  the  temperature  and 
X  is  a  variable  sampled  from  a  Chi-Square  distribution  with  v 
degrees  of  freedom.  (A  separate  sampling  is  made,  of  course, 
for  each  molecule  whose  internal  energy  is  to  be  assigned.) 

This  relation  applies  to  the  initialization  of  molecules  enter¬ 
ing  the  solution  region  through  the  sonic  orifice.  The  purpose 
of  Eq.  (2)  is  merely  to  illustrate  that  the  sampling  described 
in  the  last  yearly  report  for  initial  internal  energy  levels 
was  essentially  a  sampling  from  a  Chi-Square  distribution, 
although  it  was  not  labeled  as  such. 


3.3  Inelastic  Bimolecular  Collisions 

The  relations  used  to  define  the  post-collision  state 
vectors  following  an  inelastic  collision  (Section  4.2  of  Ref. 

1)  involve  summing  the  total  energy  of  the  collision  (internal 
energy  of  the  two  molecules  plus  the  translational  energy  of 
their  relative  motion)  and  then  redistributing  it  according  to 
the  number  of  degrees  of  freedom  for  each  of  the  three  contribu¬ 
ting  parts.  The  procedure  involves  two  samplings  from  a  distri¬ 
bution  function,  S,  of  the  form: 


SUX) 


1 

B(v1/2,v2/2) 


(v./2-l)  ( v~/2-l) 

5  1  ( l-K x) 


(3) 


where  B  denotes  the  beta  function.  Equation  (3)  statistically 
describes  the  partitioning  of  a  given  amount  of  energy  between 
two  modes,  where  v.  and  are  the  number  of  degrees  of  freedom 


of  the  two  modes,  and  the  fraction  of  the  energy  going  to 

the  mode  with  degrees  of  freedom.  Equation  (3)  was  first 
sampled  to  partition  the  total  energy  between  the  translational 
and  internal  modes,  and  then  sampled  again  to  partition  the  total 
internal  energy  between  the  two  molecules. 

The  sampling  of  Eq.  (3)  was  previously  done  via  the 
acceptance-rejection  method,  and  this  was  complicated  by  the 
presence  of  two  parameters  (v^  and  V2) ,  requiring  that  normali¬ 
zation  constants  either  be  stored  in  two  dimensional  arrays  or 
that  they  be  recomputed  at  every  sampling. 

As  discussed  in  Section  4.2,  the  sampling  of  Eq.  (3)  can 
be  achieved  in  a  more  basic  fashion  via 

X1 

^1  +  X2  ' 

where  X.^  is  selected  from  a  Chi-Square  distribution  with  v.^ 
degrees  of  freedom  and  X2  is  selected  from  a  Chi-Square  distri¬ 
bution  with  v2  degrees  of  freedom.  This  relation  uncouples  the 
sampling,  so  that  it  reduces  the  case  of  sampling  from  a  two 
parameter  distribution  to  that  of  sampling  twice  from  a  one 
parameter  distribution.  The  latter  approach  is  generally  to 
be  preferred,  since  optimization  of  the  sampling  routine  is 
much  easier  for  only  one  variable. 

Equation  (4)  is  a  lot  more  signficant  than  merely  indicating 
a  new  way  to  sample  from  the  necessary  distribution.  As  was 
shown  in  the  previous  section,  a  Chi-Square  distribution  describes 
the  energy  allocation  that  is  to  be  expected  in  an  equilibrium 
gas.  Hence,  Eq.  (4)  indicates  that  it  is  proper  to  first  select 
Chi-Square  values  as  would  be  done  for  an  equilibrum  gas,  and 
then  normalize  the  Chi-Square  values  selected  for  the  various 
competing  modes  so  that  the  available  energy  in  a  collision  is 
precisely  conserved. 


This  concept  leads  directly  to  an  even  simpler  method  for 
determining  the  post-collision  energy  distribution,  namely  each 
of  the  three  competing  modes  is  given  a  fraction  of  the  total 
energy,  given  by 

Xi 

^i  =  Xx  +  X2  +  X3  '  (5) 

where  three  Chi-Square  samplings,  each  from  a  distribution  with 
the  appropriate  number  of  degrees  of  freedom,  replace  the  two 
samplings  of  Eq.  (3) .  It  is  a  relatively  simple  matter  to  show 
that  Eq.  (5) ,  arrived  at  here  intuitively,  is  in  fact  formally 
correct. 

3.4  Equilibrium  Collision  Aftermath 

One  of  the  principal  historic  drawbacks  of  the  direct  simu¬ 
lation  Monte  Carlo  method  is  its  inefficiency  (as  opposed  to 
invalidity)  as  collisions  become  more  dominant  and  the  flow 
approaches  equilibrium.  The  necessity  of  sampling  a  very  large 
number  of  collisions  is  time  consuming  and  somewhat  redundant 
since,  once  equilibrium  is  achieved,  further  collisions  have  no 
effect  on  the  velocity  distribution.  This  concept  has  been 
utilized  in  the  past  to  cut  off  the  sampling  of  collisions  after 
the  simulation  of  a  sufficient  number  to  guarantee  equilibrium 
(Ref.  2) ,  but  it  still  required  the  simulation  of  a  large  number 
of  collisions  and  therefore  resulted  in  a  relatively  inefficient 
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simulation.  Since  the  equilibrium  limit  is  precisely  the  limit 
in  which  continuum  fluid  mechanic  descriptions  become  valid, 
Monte  Carlo  simulations  have  not  been  extensively  used  in  highly 
collisional  situations. 


Using  the  relations  of  this  section,  with  some  extensions, 
it  is  possible  to  bypass  collision  sampling  altogether  ior  cells 
in  equilibrium.  New  velocity  and  internal  energy  elements  to 
the  state  vectors  can  be  selected  directly  from  the  distributions 
resulting  from  many  collisions  (i.e.,  equilibrium)  with  the  con¬ 
straints  that  total  momentum  and  energy  be  precisely  conserved. 
The  result  is  that  cells  in  equilibrium  (i.e.,  those  near  the 
orifice  in  a  mass  spectrometer)  can  potentially  become  the 
easiest  cells  to  simulate  rather  than  the  hardest;  and  the  over¬ 
all  simulation  efficiency  can  be  substantially  improved. 

The  procedure  for  sampling  from  the  equilibrium  aftermath 
of  many  collisions  is  as  follows: 

1)  Evaluate  the  following  sums  over  all  molecules  in 
the  cell: 


-  £ 


Vi 


•  £ 


W.m.u. 

111 


=  £ 


W.m. v. 
111 


=  £ 


Wmw. 

111 


£  W.m. (u,2  +  v,2  +  w,2) 


(10) 


where  VT  ,  itk  and  E^  are  the  statistical  weighting 
factor,  the  mass  and  the  internal  energy,  respectively, 
of  the  ith  molecule,  and  u^,  and  wi  are  its  velocity 
components . 

Compute  the  center  of  mass  velocity  components,  u*,  v* 


and  w* 

u*  = 

via: 

VS1 

9 

(12) 

v*  = 

VS1 

9 

(13) 

and 

w*  = 

VS1 

• 

(14) 

The  total  translational  energy  of  the  relative  motion 
between  the  molecules,  Etrn,  can  be  represented 


Etrn  =  I  2  Wimi[(ui“u*)2+(vi-v*)2+(wi-w*)2]  ,  (15) 


although  it  is  more  easily  evaluated  via  the  equiva¬ 
lent  expression 


(16) 


The  total  cell  energy,  E.  . , 


is  therefore  given  by 


E  u*,  v*  and  w*  are  the  quantities  which  are  to  be  conserved 
in  the  equilibrium  sampling  of  new  velocity  and  internal  energy 
values  for  the  cell  molecules.  The  basic  concept  is  to  sample  a 
Chi-Square  value  appropriate  to  each  available  energy  mode,  and 
then  allocate  the  available  energy  in  proportion  to  the  assigned 
Chi-Square  values.  Sampling  a  Chi-Square  value  for  the  internal 
energy  of  each  molecule  is  straightforward,  but  the  translational 
mode  requires  a  little  care.  Since  it  is  the  energy  of  relative 
motion  that  we  are  considering,  and  a  molecule  can  have  no  rela¬ 
tive  motion  with  respect  to  itself,  all  but  one  of  the  molecules 
should  have  a  Chi-Square  value  sampled  for  its  translational 
degrees  of  freedom.  (These  Chi-Square  values  are  sampled  for 
the  three  translational  degrees  of  freedom.)  Velocity  compon¬ 
ents  of  the  molecules  will  be  assigned  one  molecule  at  a  time, 
and  when  this  is  done  the  center  of  mass  velocity  components  of 
the  remaining  molecules  are  implied  via  conservation  of  u* ,  v*, 
and  w*.  When  there  is  only  one  remaining  molecule  this  means 
that  its  velocity  components  are  implied  by  the  choices  made 
for  the  others,  and  it  does  not  have  independent  translational 
degrees  of  freedom.  The  identity  of  the  "last"  molecule  is 
arbitrary  and  has  no  effect  on  the  assigned  velocity  components. 

Let  represent  the  Chi-Square  value  sampled  for  the  ith 
molecule's  internal  energy  mode,  and  Xt^  the  value  sampled  for 
its  translational  mode.  (The  last  molecule  is  assigned  Xfci  =  0.) 
The  equilibrium  post-collision  sampling  continues  as  follows: 

5)  A  weighted  sum  of  Chi-Square  values  is  defined  via: 

s7  =  2  +  *ti)  •  (18) 

6)  The  first  molecule  is  assigned  an  internal  energy 
given  by 


(19) 


P  _  II  tot 

EI1  ■  — —  ' 

and  a  translational  energy  given  by 

„  _  XtlEtot 

Eti  "  — si —  * 


(20) 


The  relative  speed,  q^,  between  the  first  molecule  and 
the  center  of  mass  velocity  of  all  molecules  (including 
itself)  which  corresponds  to  Efcl  is  given  by 


qr  =  V2Eti(si  “  wim1)/(s1m1)  *  (21) 

The  direction  for  the  relative  velocity  is  selected  at 
random,  giving  the  three  relative  velocity  components 

url'  vrl  and  wrl  via 


A  = 

1  -  2R  , 

(22) 

B  = 

qrVl  -  A?  , 

(23) 

C  = 

2ttR  , 

(24) 

url 

=  qrA  , 

(25) 

vrl 

=  Bcos(C) 

(26) 

and 

wrl 

*  Bsin(C)  ; 

(27) 

where  each  appearance  of  R  denotes  a  distinct  evalua¬ 
tion  of  a  random  variable. 

The  updated  velocity  components  for  the  first  molecule 
are  then  given  simply  by: 
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I.\. 


U1  =  u  +  uri  . 


Vi  =  v  +  vrl 


(28) 

(29) 


and 


w, 


=  w  +  w 


(30) 


'1  rl 

10)  The  first  molecule's  contribution  to  u*,  v*  and  w*  is 
then  removed  via  the  replacements: 


<u*>new 

=  u*  -  Furl  , 

(31) 

<v*>new 

=  V*  -  Fvrl 

(32) 

and 

(w*) 

new 

-  V*  -  Fwrl  ; 

(33) 

where 

W,  m. 

F  =  _ 

1  1 

(34) 

h 

'  Wlml  * 

The  first  molecule's  contribution  to  S. 

is  then  removed 

via  the 

replacement 

(S. ) 

1  new 

=  S1  ”  wimi 

(35) 

12)  Steps  6  through  11  are  then  repeated  for  each  molecule 
in  the  cell,  except  that  when  trie  last  molecule  is 
reached  its  velocity  components  are  simply  the  cen  :er 
of  mass  velocity  components,  as  discussed  above. 


Although  the  description  of  equilibrium  sampling  with 
strict  conservation  of  total  mass,  moment  am  and  energy  may  have 
seemed  somewhat  long,  the  relations  are  quite  fast  computationally, 
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and  the  increased  efficiency  over  the  sampling  of  many  individual 
collisions  can  be  substantial.  The  degree  of  increased  effi¬ 
ciency  depends,  of  course,  on  the  collision  rate  in  a  particular 
ce.l. 

However,  some  caution  must  be  exercised  in  the  use  of 
the  equilibrium  relations.  They  are  only  justified  when 
the  number  of  collisions  simulated  in  a  given  time  step  would 
be  sufficient  to  guarantee  equilibrium.  For  the  translational 
mode  this  is  typically  on  the  order  of  three  collisions  per 
molecule,  but  since  internal  modes  (usually  rotation  for  the 
mass  spectrometer  problem)  are  also  equilibrated,  the  number  of 
collisions  should  be  somewhat  higher  (typically  five  or  so) . 
Hence,  the  major  disadvantage  in  using  the  equilibrium  relations 
is  the  decision  process  required  to  use  it  only  when  justified. 
For  that  reason,  these  relations  are  not  currently  implemented 
in  the  mass  spectrometer  code,  although  they  have  been  coded 
and  could  be  inserted  whenever  required  by  a  particularly  col¬ 
lision  dominated  case. 

4.  CHI-SQUARE  PROBABILITY  DISTRIBUTION 
4.1  Physical  Basis 

Fundamentally ,  the  Chi-Square  function  represents  the  dis¬ 
tribution  of  energy  in  an  equilibrium  classical  system  with  v 
degrees  of  freedom.  It  is  a  well  known  classical  result  that 
each  degree  of  freedom  for  a  molecule  in  an  equilibrium  gas  will 
have,  on  the  average,  an  energy  of  kT/2,  where  k  is  Boltzmann's 
constant  and  T  is  temperature.  (For  example,  the  translational 
mode,  with  three  degrees  of  freedom,  has  an  average  energy  of 
3kT/2  per  molecule.  The  distribution  of  translational  energy 
among  the  various  molecules  follows  a  Chi-Square  distribution 


with  3  degrees  of  freedom.)  Other  modes  of  energy  (molecular 
rotation  and  vibration)  have  their  own  characteristic  number  of 
degrees  of  freedom,  which  may  or  may  not  be  fully  excited  in  the 
energy  range  of  interest.  If  a  mode  is  not  fully  excited,  that 
simply  means  that  it  is  behaving  as  if  it  had  a  non-integer 
number  of  degrees  of  freedom,  within  the  classical  approximation . 
The  number  of  internal  degrees  of  freedom  is  directly  related  to 
the  heat  capacity  of  the  gas  and,  essentially,  v  is  selected  to 
match  the  known  heat  capacity  of  a  given  molecule  in  a  given 
energy  range.  The  assumption  of  a  constant  number  of  degrees 
of  freedom  is  therefore  equivalent  to  the  assumption  of  a  con¬ 
stant  heat  capacity.  A  discussion  of  the  implementation  of  such 
a  model  allowing  for  a  finite  rate  relaxation  towards  equilibrium 
between  translational  and  internal  modes  is  given  in  Ref.  1. 

4.2  Definition  and  Mathematical  Properties3 

The  Chi-Square  probability  density  function,  f(X;v),  defines 
a  distribution  of  X  in  a  domain  of  zero  to  infinity  via 

X ( V/2  “  1)exp(-X/2)  , 


where  v  is  a  positive  parameter  of  the  distribution  referred  to 
as  the  number  of  degrees  of  freedom.  The  Chi-Square  distribution 
results  in  a  mean  value  of  X  equal  to  v.  Figure  1  is  a  plot  of 
the  Chi-Square  probability  density  function  for  v  equal  to  1,  2 
and  3. 


Abramowitz,  M. ,  and  Stegun,  I. A.,  Handbook  of  Mathematical 
Functions,  National  Bureau  of  Standards,  1968,  pp.  940,  944. 


The  Chi-Square  distribution  has  a  fundamental  addition 
property  such  that  if  X^  is  selected  from  a  Chi-Square  distribu 
tion  with  degrees  of  freedom,  and  X2  is  selected  from  a  Chi- 
Square  distribution  with  v2  degrees  of  freedom;  then  their  sum 
will  be  distributed  according  to  a  Chi-Square  distribution  with 
+  v2  degrees  of  freedom.  This  property  is  of  substantial 
theoretical  and  practical  importance. 

If  the  variable  Z  is  distributed  according  to  a  normal 
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distribution  with  zero  mean  and  unit  variance,  then  Z  will 
follow  a  Chi-Square  distribution  with  one  degree  of  freedom. 

It  follows  from  the  above  addition  property  that,  in  general, 
if  Z^,  Z2,  ...,  Zn  are  n  variables  selected  from  such  a  normal 
distribution,  and  X  is  defined  as  the  sum  of  the  squares  of  the 
Z^,  then  the  X's  that  result  will  be  distributed  according  to  a 
Chi-Square  distribution  with  n  degrees  of  freedom. 

Finally,  if  t  is  distributed  according  to  a  probability 
desnity  function  g(t;p,q),  where 


g(t;p,q) 


tP'^l-t)^1 

B(p,q) 


B (p,q) 


X 

=  J  tp-1(l-t)q-1 


at  r(p+q) 


(B  is  the  Beta  function)  then  t  can  be  sampled  via 


t  = 


X1  +  X2 


where  X^  is  selected  from  a  Chi-Square  distribution  with 
degrees  of  freedom,  and  X2  is  selected  from  a  Chi-Square  distri 
bution  with  v->  degrees  of  freedom,  with 


The  significance  of  Eq.  (39)  is  that  it  reduces  the  sampling  from 
a  two  parameter  distribution  (Eq.  (37))  to  two  samplings  from  a 
one  parameter  distribution.  The  distribution  represented  by 
Eq.  (37)  arises  in  cases  where  a  constrained  amount  of  total 
energy  is  distributed  among  various  modes,  and  its  relation  to 
the  Chi-Square  distribution  apparently  has  not  been  appreciated 
by  developers  of  techniques  for  Monte  Carlo  fluid  mechanics. 

4.3  Sampling  From  a  Chi-Square  Distribution 

The  need  for  sampling  from  a  Chi-Square  distribution  comes 
up  when  sampling  initial  values  of  internal  energies,  when 
calculating  inelastic  collisions  via  the  statistical  collision 
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model  or  when  calculating  the  equilibrium  aftermath  of  many 
collisions  in  a  cell.  Since  these  operations  must  be  performed 
repeatedly  in  the  heart  of  a  Monte  Carlo  simulation,  it  is  import 
ant  that  the  sampling  be  done  efficiently  and  accurately. 

For  clarity,  the  result  of  each  sampling  method  discussed 
below  will  be  denoted  by  a  different  letter  subscript  to  X. 

All  sampling  procedures  make  use  of  a  random  number  generator 
which  returns  a  number,  R,  selected  from  a  probability  density 
which  is  uniform  on  the  interval  between  zero  and  one.  Each 
occurrence  of  R  indicates  a  distinct  sampling  from  the  random 
number  generator. 

4 

Borgnakke,  Claus,  and  Larsen,  Paul  S.,  "Statistical  Collision 
Model  for  Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture," 
Journal  of  Computational  Physics,  Vol.  18,  1975,  pp.  405-420. 


4.3.1  Analytic  Sampling  for  Integer  v 

Direct  sampling  of  Eq.  (36)  can  be  performed  for  integer 
v,  as  shown  below. 

4. 3. 1.1  v  =  0 

As  v  (an  intrinsically  nonnegative  quantity)  approaches 
zero,  the  distribution  function  approaches  a  delta  function, 
and  a  proper  sampling  is  achieved  by  simply  selecting 

Xa  =  0  .  (42) 


4. 3.1.2  v  =  1 


For  sampling  with  v  =  1  (as  well  as  for  several  other  cases) 

•  .  .  2 

it  is  convenient  to  introduce  the  transformation  Z  =  X.  Z  is 

then  distributed  according  to  the  probability  density  function 

p(Z)  given  by 


P(Z) 


Z (v~1)exp(-Z2/2) 
2(v/2  -  1)r(v/2) 


(43) 


For  v  —  1 ,  this  distribution  is  simply  a  normal  distribution 
adjusted  to  allow  for  positive  only  argument.  Sampling  from 
this  distribution  is  described  in  Ref.  3.  When  the  result  is 
cast  back  in  terms  of  X,  the  result  is 


A  =  2ttR 

and 


(44) 


X,  =  -21og(R)  sin2  (A) 


(45) 


4. 3. 1.3  v  =  2 


When  v  =  2,  the  integral  of  Eq.  (36)  can  be  analytically 
inverted,  leading  to  the  direct  sampling 

Xc  =  -21og (R)  .  (46) 


4. 3. 1.4  v  Equal  to  an  Even  Integer 

The  extreme  simplicity  of  the  above  sampling  for  v  =  2, 
together  with  the  addition  property  of  the  Chi-Square  distribu¬ 
tion,  means  that  sampling  for  v  equal  to  an  even  integer  is  quite 
direct.  Let  J  =  v/2,  then  a  proper  Chi-Square  sampling  is  given 
by 

Xd  =  -21og(R1R2...Rj)  ,  (47) 

where  R^  through  Rj  denote  j  samplings  from  the  random  number 
generator.  The  fact  that  the  log  need  only  be  taken  once  in 
Eq.  (47)  means  that  the  evaluation  of  Xd  is  quite  efficient, 
even  for  moderately  large  v. 


4. 3. 1.5  v  Equal  to  an  Odd  Integer 

For  v  equal  to  an  odd  integer,  the  addition  property  of  the 
Chi-Square  distribution  allows  the  simple  combination  of  the 
results  for  v  equal  to  one  and  v  equal  to  an  even  integer,  i.e., 

Xe  =  Xb  +  Xd  '  (48) 

where  X^  is  given  in  Eq.  (45)  and  Xd  is  given  in  Eq.  (47)  with 
J  =  (v-D/2. 
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4.3.2  Generalized  Acceptance-Rejection  Sampling 


For  non-integer  v,  it  is  necessary  to  use  a  generalized 
form  of  acceptance-rejection  sampling.  Before  the  application 
to  Chi-Square  sampling  is  presented,  the  acceptance-rejection 
technique  and  its  generalization  will  be  briefly  discussed. 


4. 3. 2.1  Standard  Acceptance-Rejection  Sampling 


The  usual  acceptance-rejection  technique  for  sampling  from 
a  general  distribution  function,  p(x),  proceeds  as  follows: 


1)  The  domain  of  x  is  approximated,  if  necessary,  by  a 
finite  sub-domain. 


2)  The  maximum  value  of  p(x),  p*,  is  calculated. 


3)  A  variable  £  is  selected  from  the  domain  of  x  via 


^  xmin  +  R^xmax  ~  xmin^ 


4)  p(U/p*  is  calculated,  and  another  random  variable, 
R  is  generateo.  x  is  set  equal  to  £  if  R  is  less 
than  p(£)/p*. 


5)  Steps  3  and  4  are  repeated  until  a  value  of  x  is 


determined. 


Note  that  the  probability  of  acceptance  of  the  random  variable 
in  step  4  is  proportional  to  the  distribution  function  being 
sampled,  so  the  resulting  x  values  will  follow  the  desired 
distribution  function. 


Although  the  generality  of  this  approach  makes  it  very 
powerful,  it  does  suffer  from  the  following  drawbacks: 


If  the  distribution  function  differs  significantly 
from  its  maximum  value  within  a  substantial  portion  of 
the  sampled  domain,  then  the  rejection  rate  may  be 
high.  Thic  obviously  .cuds  to  a  al-j»  s»ainplj.na  procedure 
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•  If  the  finite  sub-domain  is  reduced  to  increase  the 
acceptance  rate,  then  the  sampling  deviates  from  the 
true  distribution  function. 

•  The  procedure  is  incapable  of  sampling  from  an 
unbounded  distribution  function. 

4. 3. 2. 2  Generalization  of  the  Acceptance-Rejection  Technique 

The  following  procedure  comprises  a  generalization  of  the 
acceptance-rejection  technique: 

1)  A  second  distribution  function,  q(x) ,  which  can  be 
sampled  analytically  is  chosen.  Conditions  on  q(x) 
will  be  discussed  below. 

2)  The  maximum  value  of  p(x)/q(x),  (p/q)*,  is  calculated. 

3)  A  variable,  £,  is  sampled  from  q. 

4)  Q  =  [p(£)/q(£)]  /(p/q)  *  is  calculated,  and  another 
random  variable,  R,  is  generated,  x  is  set  equal 
to  E,  if  R  is  less  than  Q. 

5)  Steps  3  and  4  are  repeated  until  a  value  of  x  is 
determined . 

It  should  be  noted  that  the  probability  density  for  a  given 
value  of  x  is  proportional  to  the  product  of  the  initial  selec¬ 
tion  probability  times  the  acceptance  probability.  Since  the 
former  probability  is  proportional  to  q(x),  and  the  latter  is 
proportional  to  p(x)/q(x),  the  distribution  of  accepted  values 
does  indeed  follow  the  distribution  function  p(x). 

The  usual  acceptance-rejection  technique  is  simply  the 
case  where  q(x)  is  constant,  but  it  is  evident  that  this  is  not 
always  the  best  (or  even  a  possible)  choice.  All  of  the 
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objections  to  the  standard  acceptance-rejection  technique  can 
be  removed  or  ameliorated  by  a  suitable  choice  for  q(x).  In 
particular : 

•  There  is  no  need  to  approximate  the  domain  of  x  with 
a  finite  sub-domain.  It  is  merely  necessary  that  the 
domain  for  q  include  the  domain  for  p.  The  domain 
for  q  can  be  larger  than  that  for  p,  since  whenever  a 
value  is  selected  from  outside  the  domain  for  p  it 
will  always  be  rejected  in  step  4  above. 

•  If  q  is  selected  to  be  close  to  p,  at  least  in  the 
region  of  highest  probability,  then  the  acceptance 
rate  of  trial  values  will  be  large. 

•  Unbounded  distribution  functions  can  be  sampled  if 
q  is  chosen  to  have  the  same  type  of  singularity  as 
p,  since  the  only  requirement  is  that  the  ratio  (p/q) 
remain  bounded. 

For  any  given  situation,  the  choice  of  the  function  q  is 
a  bit  of  an  art,  guided  by  the  concerns  highlighted  above:  q 
must  have  a  domain  which  includes  the  domain  of  p;  p/q  must 
remain  bounded;  and  (p/q)  should  achieve  its  maximum  in  the 
vicinity  of  the  maximum  of  p. 

4.3.3  Exact  Acceptance-Rejection  Sampling  for  a  Chi-Square 
Distribution  with  Large  v 

The  acceptance-rejection  technique  described  above  can  be 
used  to  achieve  an  exact  sampling  from  a  Chi-Square  distribution 
for  large  v.  (Actually,  the  approach  is  perfectly  valid  for 
all  v  >  1,  but  the  method  to  be  described  in  Section  4.3.5  is 
to  be  preferred  for  v  <  45,  or  so.)  The  procedure  utilizes  the 
transformed  Chi-Square  distribution,  p(Z) ,  given  by  Eq.  (43)  as 


the  distribution  to  be  sampled.  A  normal  distribution  is  used 
as  the  initial  distribution  which  can  be  sampled  analytically. 
The  normal  distribution  is  chosen  to  have  a  unit  variance  and  a 
mean  which  corresponds  to  the  location  of  the  maximum  of  p(Z). 
This  maximum  occurs  at  Z*  given  by 

Z*  =  y  v  -  1  .  (49) 

The  functional  form  of  the  normal  distribution,  q(Z),  is 

q  (Z)  =  exp[-(Z  -  Z*)2J/VTtT  ,  (50) 

which  not  only  has  a  maximum  at  the  same  location  as  Eq.  (43) , 
but  has  the  same  exponential  factor  as  Z  approaches  infinity  and 
a  domain  which  includes  that  of  Eq.  (43) .  The  sampling  of  a 
Chi-Square  value  proceeds  as  follows: 

1)  z*  =  "Vv  -  1  is  calculated. 

2)  A  sample  from  the  distribution  given  by  Eq.  (51) 


is  taken  via: 

A  =  2ttR  (51) 

B  =  -log(R)  (52) 

Z  =  Z*  +  V2Bsin(A)  (53) 

3)  The  acceptance  probability,  Q,  is  computed  as 

Q  =  (Z/Z*)  (V  “  1)exp[-Z*(Z  -  Z*)]  .  (54) 

(Q  is  taken  to  be  zero  for  negative  Z.) 


4)  Another  random  variable  is  generated,  and  Z  is  kept 

if  Q  >  R.  If  Z  is  rejected,  then  steps  2-4  are  repeated 
until  a  Z  value  is  accepted. 


5)  When  a  Z  value  is  accepted,  then  the  corresponding 
Chi-Square  value  is  given  by 

Xe  =  Z2  .  (55) 

This  procedure  is  illustrated  in  Fig.  2  which  shows  p(Z),  q(Z) 
and  Q ( Z )  for  v  =  50.  Note  that  the  acceptance  probability  is 
near  unity  in  the  vicinity  of  the  maxima  of  the  two  distribution 
functions,  so  a  large  fraction  of  the  selected  samples  of  q(Z) 
will  be  accepted  as  samples  of  p(Z). 


4.3.4  Exact  Acceptance-Rejection  Sampling  for  a  Chi-Square 
Distribution  with  (0  <  v  <  2) 

For  this  domain  of  v  it  is  convenient  to  introduce  another 
transformation  to  Eq.  (36) .  If  W  is  defined  by 

W  =  exp (-X/2)  ,  (56) 

then  the  probability  density  function  for  W  is  given  by  h(W) , 
where 


h(w)  = 


[-log(W)]<v' 

r(v/2) 


(v/2  -  1) 


The  domain  for  W  is  finite  (between  0  and  1),  but  h(W)  becomes 
infinite  as  W  approaches  unity.  The  generalized  acceptance- 
rejection  technique  can  still  be  used,  however,  since  the  func¬ 
tion  q(W)  given  by 


q(W)  =  (v/2) (1  -  W) {V/2  ”  1} 


(58) 


has  the  same  type  of  singularity  and  can  be  analytically  sampled 
The  Chi-Square  sampling  for  (0  <  v  <  2)  proceeds  as  follows: 


Figure  2.  A  Representation  of  the  Transformed  Chi-Square 
Distribution,  p(Z),  for  v  =  50  (solid  line).  p(Z)  is 
sampled  by  first  selecting  a  variable  from  the  shifted 
normal  distribution  (dashed  line)  and  keeping  it  with  a 
probability  given  by  Q(Z)  (dotted  line) . 
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1)  A  sample  from  q(W)  is  generated  via 

W  =  1  -  R(2/v)  .  (59) 

2)  The  acceptance  probability,  Q,  is  computed  from 

Q  =  [(W  -  l)/log(W)]  {1  "  v/2)  .  (60) 

3)  Another  random  variable,  R,  is  generated,  and  W  is 
kept  if  Q  >  R;  otherwise  steps  2  and  3  are  repeated 
until  a  value  for  W  is  accepted. 

4)  When  a  value  for  W  is  accepted,  the  corresponding 
Chi-Square  value  is  given  by 

Xf  =  -21og (W)  .  (61) 

This  procedure  is  illustrated  in  Fig.  3,  which  shows  the  two 
distribution  functions,  h(W)  and  q(W),  and  the  acceptance  prob¬ 
ability,  Q(W) ,  for  v  =  1.  It  can  be  seen  that  q(W)  provides  an 
excellent  choice  for  the  initial  selection  of  W,  since  the 
acceptance  probability  remains  high  throughout  the  important 
domain  of  W.  This  point  will  be  discussed  in  more  detail  in 
Section  4.3.6. 

4.3.5  Exact  Chi-Square  Sampling  for  General  v 

Using  the  fundamental  addition  property  of  Chi-Square 
distributions,  it  is  possible  to  combine  the  procedure  described 
in  Section  4. 3. 1.4  for  v  equal  to  an  even  integer  with  the  pro¬ 
cedure  described  in  Section  4.3.4  for  (0  <  v  <  2)  to  achieve  an 
exact  general  sampling  technique  for  arbitrary  v.  This  is  given 
simply  by 

Xg  =  Xd  +  Xf  ,  (62) 
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q(W)  (Eq.  58) 
RCCEPTRNCE  PROBRBILITY 


where  X^  is  calculated  from  Eq.  (47)  with  J  equal  to  the  integer 
portion  of  v/2,  and  is  calculated  as  in  the  preceding  section 
with  v  being  replaced  by  v  -  2J. 

It  is  to  be  noted  that  both  the  approach  given  in  this 
section  (Eq.  (63))  and  that  given  in  Section  4.3.3  (Eq.  (55)/ 
are  exact  and  applicable  for  v  >  1.  In  general,  the  approach  of 
this  section  is  considerably  faster,  although  as  v  gets  large 
the  approach  of  Section  4.3.3  becomes  more  attractive.  There 
are  two  potential  difficulties  with  Eq.  (62)  as  v  becomes  very 
large.  Firstly,  the  product  required  in  Eq.  (47)  gets  more  and 
more  cumbersome  to  compute  as  v  increases,  and,  secondly,  the 
larger  the  number  of  factors  in  this  product  the  greater  is  the 
chance  that  it  will  yield  a  number  so  small  as  to  produce  a 
floating  point  underflow  on  a  computer.  (Since  Monte  Carlo 
codes  must  be  highly  reliable,  any  such  problem  should  be  made 
essentially  impossible.)  It  turns  out  that  the  second  problem 
is  more  restrictive  (at  least  for  32  bit  computers) ,  dictating 
that  the  Eq.  (55)  should  be  used  for  v  greater  than  45,  or  so. 
This  keeps  the  probabilityof  an  underflow  below  10  10  on  any 
given  sampling. 

4.3.6  Approximate  Chi-Square  Sampling  for  (0  <  v  <  2) 

The  procedure  described  in  the  preceding  section  is  quite 
efficient,  but  it  is  nonetheless  useful  to  consider  approximate 
methods  for  sampling  from  Ch.> -Square  distributions.  While  it 
would  be  scarcely  possible  to  improve  on  sampling  for  even 
integer  v  discussed  in  Section  4. 3. 1.4,  it  is  reasonable  to 
investigate  approximations  for  the  (0  <  v  <  2)  portion  discussed 
in  Section  4.3.4.  A  likely  place  to  look  for  useful  approxima¬ 
tions  in  this  procedure  is  in  the  calculation  of  Q  (Eq.  (60)), 


which  must  be  performed  for  every  W  selected  in  Eq.  (59) .  (Note 
that  the  calculation  of  Q  involves  more  computational  effort 
than  the  calculation  of  W.) 

The  overall  probability  that  the  value  chosen  in  Eq.  (59) 
will  be  kept  as  a  sample  of  Eq.  (57)  is  given  by  P,  where 
1 

P  =  f  q(W)Q(W)dW  =  r (1  +  v/2)  .  (63) 

0 

Hence,  as  v  approaches  0  or  2 ,  all  initially  selected  values  of 
W  are  kept  as  valid  samples  of  Eq.  (57)  and  the  computation  of 
Q  serves  no  useful  purpose.  In  the  worst  case  (v  =  0.92)  the 
overall  acceptance  probability  is  89%,  and  only  11%  of  the 
initially  selected  variables  are  rejected.  The  approximate 
Chi-Square  sampling  involves  approximating  0  (W)  by  an  easily 
calculable  function  which  differs  little  from  Eq.  (60).  The 
current  approximation  is  <0  ('.*.’)  ,  riven  bv 

a 

Qa(W)  =  1  -  (1  -  ./2M1  -  -i)[.5  -  ..Mi  -  W)2]  ,  (64) 

where 

a  (v)  =  .2511v  +  .2073  .  ,65) 

Q  was  selected  to  match  the  value  and  slope  of  Q  at  W  =  1, 
which  is  the  region  of  highest  probability  density.  The  coeffi¬ 
cient  a(v)  is  a  linear  fit  to  values  chosen  to  be  optimal  in  the 
least  squares  sense.  A  comparison  of  Q  and  Q  for  v  =  1.0  is 
shown  in  Fig.  4,  which  demonstrates  the  substantial  accuracy  of 
the  approximation. 

It  is  fundamentally  more  important,  of  course,  to  compare 
the  correct  Chi-Square  distribution  with  the  distribution  which 
is  effectively  being  sampled  in  the  approximate  technique.  If 


h3 (W)  is  the  approximation  analog  of  h(W),  then  h  (W)  is  pro- 
a  a 

portional  to  the  product  q(W)Q(W),  i.e., 

ha(W)  «  A (1  -  W)  (v/2  "  1}Qa(W)  ,  (66) 

where  the  normalization  factor.  A,  is  determined  by  requiring 
that  ha (W)  give  unity  when  integrated  between  zero  and  one. 

This  results  in 


1 

A 


2 

v  +  2v  +  8 
2v2  +  4v 


2a  (v) 
v  +  6 


(67) 


Once  h  (W)  is  defined,  the  corresponding  distribution  function 

CL 

for  X,  f  (X; v)  is  obtained  by  multiplying  h  (W)  by  the  magnitude 
of  dW/dX  (-VI/2)  ,  and  substituting  W  =  exp(-X/2).  The  comparison 
between  f(X;1.0)  and  f  (X;1.0)  is  given  in  Fig.  5,  and  the 

Cl 

agreement  is  excellent.  The  use  of  the  approximate  technique  is 
approximately  40%  faster  than  the  exact  acceptance-rejection 
technique,  and  the  difference  in  the  distributions  being 
sampled  will  probably  always  be  negligible.  Although  the 
ability  to  sample  from  an  exact  Chi-Square  distribution  will  be 
kept  as  an  option,  it  is  felt  that  the  approximate  technique 
offers  a  substantial  time  savings  for  an  inconsequential  loss 
of  accuracy. 


5.  INCLUSION  OF  ELECTRIC  FIELD  EFFECTS 
5.1  Acceleration  of  Charged  Species 

An  essential  element  of  the  computational  model  is  the 
accelerating  effects  of  electric  fields  on  charged  species. 

The  program  was  generalized  this  year  so  that  axial  position 
and  velocity  elements  of  a  molecular  state  vector  are  updated 
to  include  the  effect  of  electric  fields  in  the  portion  of  the 
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program  that  advances  the  molecules  along  their  trajectories. 
(The  term  molecule  is  used  here  in  a  general  sense,  and  is  meant 
to  apply  to  cluster  ions  as  well.)  If  a  molecule  has  an  initial 
axial  position  and  velocity  given  by  zQ  and  vzQ,  then  after  a 
time  t,  its  axial  position  and  velocity  will  be  given  by  z^  and 
v2i,  where 


=  2 , 


+  t(V 


zO 


2Et) 
2M  ' 


and 


(68) 


zl 


v  + 
zO  M 


(69) 


In  these  relations,  which  replace  Eqs.  (58)  and  (61)  of  Ref. 

1,  q  is  the  molecular  charge,  M  is  the  molecular  mass  and  E  is 
the  local  electric  field  strength.  The  electric  field  strength 
is  calculated  from 

E  =  -V4>  ,  (70) 

where  <J>  denotes  the  electric  potential.  The  electric  field  is 
assumed  to  be  one  dimensional  in  the  axial  direction,  and  is 
calculated  from  input  grid  voltages  and  positions.  (The  assump¬ 
tion  of  one  dimensionality  is  consistent  with  the  geometry  of 
the  mass  spectrometer  and  the  primary  interest  of  flow  along 
the  axis  of  symmetry.)  Note  that  for  q  =  0,  the  above  relations 
reduce  to  Eq.  (58)  and  Eq.  (61)  of  Ref.  1  as  applied  to  uncharged 
species. 


5.2  Neglect  of  Space  Charge  Contribution  to  E 


The  calculation  of  E  from  the  input  grid  conditions,  neglect 
ing  the  contribution  from  the  ions  in  the  flow,  is  a  considerable 
simplification  which  results  in  an  uncoupling  of  the  electric  and 


flow  field  problems.  The  validity  of  this  approximation  for 
the  case  of  interest  can  be  demonstrated  via  a  simplified  example 
Consider  two  grids  of  potentials  ^  and  (j^,  with  the  former  being 
located  at  2  =  0  and  the  latter  being  located  at  z  =  s.  If 
there  is  a  uniform  space  charge  density,  p,  between  the  plates, 
then  the  potential  between  the  plates  is  determined  as  the  solu¬ 
tion  to  Poisson's  equation: 


V2<f> 


d 

dz2 


(71) 


which  is 

4>  = 


$2  +  (<J>2 


+  ~  z(s 

is 


-  z) 


(72) 


(In  these  relations,  cn  is  the  permittivity  of  free  space. 


*”12 


equal  to  8.855x10  J"“  farad/m.) 
strength  is  then  computed  as 


The  corresponding  electric  field 


=  _  = 
dz 


^1  ~  ^2  2z 

- -  1  +  H(—  -  1) 

s  s 


(73) 


where  the  dimensionless  parameter,  H,  is  defined  by 

2 


H  = 


ps 


2Eq  (4>2  "  <t>2) 


(74) 


The  purpose  of  this  exercise  has  been  to  determine  the  dimension¬ 
less  parameter  which  determines  when  space  charge  has  a  signifi¬ 
cant  effect  on  the  electric  field  distribution.  That  parameter 
is  H,  and  it  can  be  seen  from  Eq.  (72)  that  the  electric  field 
distribution  is  unaffected  by  space  charge  as  long  as  H  is  small 
compared  to  unity. 


A  conservative  estimate  of  H  for  the  problem  of  interest  is 
obtained  by  taking  a  typical  potential  difference  of  10  volts 


existing  between  plates  separated  by  a  distance  on  the  order  of 
0.01  meter.  p  can  be  grossly  overestimated  by  taking  the  ambient 
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ion  density  (~10  m  )  times  the  electronic  charge  (1.6x10 
coulomb) ,  neglecting  the  mitigating  effects  of  charge  balance 

and  the  density  reduction  in  a  vacuum  expansion.  Even  with  this 

-4 

substantial  overestimate,  H  is  on  the  order  of  10  ,  and  the 

perturbation  of  the  electric  field  due  to  space  charge  is  quite 
negligible. 

6.  SUMMARY  OF  CODE  STATUS 

With  all  the  technical  detail  presented  in  this  and  the 
previous  yearly  report,  it  is  useful  to  give  an  overview  of  the 
code  current  status.  The  program  is  (temporarily,  at  least) 
called  the  EXPANDO  code,  and  it  currently  has  36  subroutines 
comprising  over  4100  lines  of  FORTRAN.  It  is  written  in  a 
highly  modular  fashion,  so  that  new  capabilities  can  be  included 
without  major  modifications  of  the  code.  There  has  been  a  sub¬ 
stantial  effort  to  make  the  code  efficient  in  its  utilization  of 
computational  resources.  As  a  result  of  this,  all  runs  to  date 
have  been  performed  satisfactorily  on  a  micro-computer.  The 
performance  on  a  larger  computer  is  naturally  expected  to  be 
even  better. 

The  code  implements  and  expands  the  direct  simulation  Monte 
Carlo  method  for  transitional  gas  dynamics,  as  applied  to  the 
flowfield  beyond  the  orifice  in  a  mass  spectrometer.  By  directly 
simulating  the  molecular  processes  of  translations,  collisions 
and  chemical  reactions  it  retains  validity  even  though  the  forma¬ 
lism  of  continuum  fluid  mechanics  breaks  down  in  this  case. 

Some  of  the  specific  features  of  the  EXPANDO  code  are: 


•  It  can  treat  an  arbitrary  number  of  species  and  an 
arbitrary  number  of  chemical  reactions  between  them. 

•  The  species  are  allowed  to  have  internal  energy,  and 
the  code  describes  the  nonequilibrium  process  by 
which  the  internal  modes  go  out  of  equilibrium  with 
the  translational  mode. 

•  Species  are  allowed  to  interact  with  a  velocity 
dependent  gas  kinetic  cross  section.  This  is  a 
particularly  important  feature  for  the  mass  spectrom¬ 
eter  problem  since  the  static  temperature  of  the  free 
jet  expansion  changes  so  radically  as  the  flow  expands 
from  the  orifice. 

•  It  can  treat  the  flowfield  between  the  orifice  and 
the  skimmer  in  one  run,  or  it  can  break  up  the  flow- 
field  into  multiple  segments  and  solve  for  them 
separately. 

•  It  can  solve  for  the  gas  dynamic  flowfield  of  the 
major  neutral  species  separately  and  then  go  back 
and  calculate  the  minor  species  solution  as  a 
perturbation. 

•  Ionic  species  are  allowed,  and  the  accelerating 
effects  of  electric  fields  produced  by  charged  grids 
can  be  simulated. 

•  The  code  performs  self  checks,  to  make  sure,  for 
instance,  that  the  collision  rate  is  being  simulated 
properly. 

Only  two  major  capabilities  must  be  added  to  the  program. 
First,  the  reflection  of  jet  molecules  off  the  skimmer  must  be 
simulated  to  make  sure  that  there  is  no  substantial  skimmer 


interference  in  the  molecular  beam.  Since  the  skimmer  will  be 
cryogenically  cooled,  the  code  will  allow  for  a  fraction  of  the 
molecules  impacting  the  skimmer  wall  to  stick.  Those  that  do 
not  stick  will  specularly  reflect  off  the  skimmer.  The  second 
new  feature  is  the  ability  to  describe  a  background  gas  within 
the  mass  spectrometer  which  may  intrude  into  the  jet  and  degrade 
the  beam  quality.  (This  background  gas  is  also  a  result  of  a 
sticking  coefficient  being  less  than  unity.) 

It  should  be  stressed  that  merely  writing  a  computer  program 
does  not  comprise  the  totality  of  the  present  effort.  There  is 
also  further  work  to  be  done  in  the  characterization  of  the 
cross  sections  for  the  critical  agglomeration  and  fragmentation 
processes  that  are  to  be  simulated.  The  computer  program  can 
only  predict  the  macroscopic  implications  of  these  cross  sec¬ 
tions;  the  generation  of  the  proper  input  requires  sound  scien¬ 
tific  judgement.  This  will  be  a  major  task  for  the  final  year 
on  this  contract. 

7.  DATA  COMPARISON  FOR  SPECIES  SEPARATION 

To  serve  as  a  test  f  the  basic  physics  and  numerics  embodied 
in  the  code,  it  was  deemed  desirable  to  make  a  calculational 
comparison  with  published  data  as  a  verification.  Naturally, 
this  comparison  can  only  involve  features  of  the  code  which  have 
already  been  instituted.  The  data  that  were  chosen  for  compari¬ 
son  involve  the  separation  of  C02  and  H2  in  a  free  jet  expansion, 
as  published  in  Ref.  5.  The  case  selected  involves  a  Reynold's 
number  which  is  in  the  relevant  range  for  the  mass  spectrometer 
flow  problem. 

5McCay,  T.  D. ,  and  Price,  L.  L. ,  "Diffusive  separation  of  binary 
mixtures  of  CO2-H2  in  a  sonic-orifice  expansion,"  Physics  of 
Fluids,  26(8),  August  1983,  pp.  2115-2119. 
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7.1  Problem  Definition 


The  case  of  interest  involves  producing  a  free  jet  expan¬ 
sion  from  a  reservoir  which  has  a  C02  mole  fraction  of  94.9%, 
with  the  remaining  5.1%  being  H2 .  The  stagnation  temperature 

was  286  K,  and  the  stagnation  pressure  is  estimated  to  be 
3  2 

1.48x10  dyne/cm  .  (It  was  necessary  to  estimate  the  stagnation 
pressure  since  this  quantity  was  not  given  by  the  experimenters. 
The  Reynold ' s  number  based  on  the  stagnation  speed  of  sound  and 
the  orifice  diameter  was  supplied,  but  the  viscosity  on  which 
the  Reynold's  number  was  based  was  not  given  either.  Since  the 
viscosities  are  well  known  for  these  gases,  the  error  induced 
here  is  probably  negligible.)  The  gas  was  expanded  through  a 
sonic  orifice  of  0.32  cm  diameter,  and  the  species  concentra¬ 
tions  were  measured  via  electron  beam  fluorescence.  The  data 
were  presented  as  the  ratio  of  CC>2  to  H2  concentrations  on  the 
jet  axis  as  a  function  of  distance  downstream  from  the  orifice. 


7.2  Computational  Considerations 

Effective  cross  sections  for  CC>2  and  H2,  together  with  their 
energy  dependence,  are  given  by  Bird.6  In  calculating  the  jet 
expansion,  C02  was  taken  to  have  2.0  internal  degrees  of  freedom, 
since  its  vibrational  energy  rapidly  freezes.  The  orifice  con¬ 
ditions  were  based  on  equilibrium,  however,  with  C02  having  3.58 
internal  degrees  of  freedom.  (H2  was  taken  to  have  2.0  internal 
degrees  of  freedom  throughout  the  calculation.)  The  discharge 
coefficient  was  calculated  using  the  relations  for  N2  as  being 


Bird,  G.  A.,  "Monte-Carlo  Simulation  in  an  Engineering  Context," 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas 
Dynamics,  Vol.  74,  Progress  in  Astronautics  and  Aeronautics, 

AIAA,  New  York,  1981. 


0.87.  (Although  this  is  certainly  not  precise,  it  is  also  an 
acceptable  approximation  for  the  present  purposes.)  The  grid 
structure  described  in  Ref.  1  was  employed,  with  the  orifice 
radius  being  divided  into  six  cells,  and  three  additional  cells 
added  beyond  the  outside  of  the  orifice  for  a  total  of  nine  cells 
in  the  radial  direction.  Subsequent  downstream  cells  were 
enlarged  to  keep  pace  with  the  expanding  jet,  until  a  total  of 
twenty  axial  stations  (180  cells  altogether)  were  utilized. 

Since  95%  of  the  flow  was  C02,  it  was  reasonable  to  perform 
this  calculation  using  the  ability  to  separately  compute  major 
and  mino_  species.  Accordingly,  the  C02  distribution  was  first 
calculated  assuming  no  H2  to  be  present.  After  the  CC>2  jet  was 
defined,  the  calculation  for  H2  diffusing  through  the  CC>2  was 
performed,  allowing  the  H2  to  suffer  collisions  with  itself  as 
well  as  C02. 

The  calculation  of  the  H2  flow  involved  some  special  prob¬ 
lems  which  were  fairly  specific  to  its  light  molecular  weight. 
Since  the  H2  is  over  20  times  lighter  than  C02,  its  thermal 
velocity  is  between  4  and  5  times  as  large  as  that  for  CC>2. 

This  has  two  implications  of  importance  for  the  calculation. 
First,  spatial  segmentation  of  this  solution  was  not  feasible. 
This  was  because  the  effective  Mach  number  of  the  H2  (considered 
by  itself)  remained  quite  small.  This  can  be  seen  by  consider¬ 
ing  an  equilibrium  mixture  flowing  at  high  speed.  If  the  flow 
is  in  equilibrium,  then  the  mean  velocity  and  temperature  of 
each  component  of  the  mixture  will  be  the  same.  For  a  lighter 
gas,  however,  a  given  temperature  translates  into  a  larger 
thermal  velocity,  and  the  Mach  number  of  the  lighter  component 
considered  alone  is  less  than  that  for  the  mixture.  The  same 
statement  is  basically  true  for  the  nonequilibrium  situation  of 
interest  here.  Although  there  were  velocity  and  temperature 
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slips  between  the  two  components  of  the  mixture,  the  effective 
Mach  number  of  the  ^  remained  small  enough  so  that  it  was  not 
reasonable  to  truncate  the  solution  region  and  neglect  backflow 
of  H2  molecules.  Therefore,  the  solution  was  carried  out  in  one 
segment  to  an  axial  distance  of  approximately  ten  orifice  diam¬ 
eters  (the  extent  of  the  data) . 

The  second  problem  of  relevance  was  the  time  step  involved. 
The  time  step  should  be  small  compared  to  the  mean  time  between 
collisions  for  a  given  molecule,  and  a  molecule  of  ,  with  its 
much  larger  thermal  velocity,  suffers  collisions  at  a  much 
faster  rate  than  a  molecule  of  CO2 •  Hence,  for  the  minor  species 
run,  the  time  step  was  reduced  by  a  factor  of  V 22 .  It  should  be 
noted  that  the  ability  to  separate  major  and  minor  species 
enabled  the  major  species  to  be  treated  with  a  much  larger  time 
step  than  the  minor  species.  If  one  solution  were  carried  out 
combining  major  and  minor  species,  then  the  entire  solution  would 
have  to  be  limited  by  the  time  step  required  by  the  H2. 

The  ability  to  compute  minor  species  with  a  different  time 
step  than  that  used  for  major  species  does  have  relevance  to  the 
mass  spectrometer  flow  problem.  In  that  problem,  the  minor 
species  of  interest  will  include  charged  species  which  are 
accelerated  by  electric  fields  relative  to  the  rest  of  the  flow. 
Such  molecules  will  have  an  elevated  collision  frequency  and 
will  therefore  require  a  smaller  computational  time  step  in 
order  to  assure  that  the  time  step  is  small  compared  to  the  mean 
time  between  collisions.  With  the  ability  to  separate  major 
and  minor  species,  the  entire  solution  need  not  be  computed 
under  this  restriction;  it  is  only  necessary  to  apply  it  when 
determining  the  minor  species  solution. 
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7.3  Results 


The  comparison  of  the  calculational  results  with  the 
published  data  is  shown  in  Fig.  6.  The  agreement  is  excellent 
This  is  felt  to  be  a  strong  confirmation  of  the  code's  ability 
to  calculate  a  basic  free  jet  flow  field. 
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Figure  6.  A  Comparison  of  Calculations  with  Published  Data 
for  the  Diffusive  Separation  of  H2  and  CO2  in  a  Free  Jet 
Expansion  From  a  Sonic  Orifice.  The  orifice  is  0.32  cm  in 
diameter,  and  the  stagnation  temperature  and  pressure  are 
286  K  and  1480  dyne/cm2. 
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